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ABSTRACT 

The structural evolution of substructure in cold dark matter (CDM) models is investigated combining 
"low-resolution" satellites from cosmological TV-body simulations of parent halos with N = 10 7 particles 
with high-resolution individual subhalos orbiting within a static host potential. We show that, as a result 
of mass loss, convergence in the central density profiles requires the initial satellites to be resolved with 
N = 10 7 particles and parsec-scale force resolution. We find that the density profiles of substructure halos 
can be well fitted with a power-law central slope that is unmodified by tidal forces even after the tidal 
stripping of over 99% of the initial mass and an exponential cutoff in the outer parts. The solution to the 
missing-satellites problem advocated by Stoehr et al. in 2002 relied on the flattening of the dark matter 
halo central density cusps by gravitational tides, enabling the observed satellites to be embedded within 
dark halos with maximum circular velocities as large as 60 kms -1 . In contrast, our results suggest that 
tidal interactions do not provide the mechanism for associating the dwarf spheroidal satellites (dSphs) 
of the Milky Way with the most massive substructure halos expected in a CDM universe. Motivated 
by the structure of our stripped satellites, we compare the predicted velocity dispersion profiles of 
Fornax and Draco to observations, assuming that they are embedded in CDM halos. We demonstrate 
that models with isotropic and tangentially anisotropic velocity distributions for the stellar component 
fit the data only if the surrounding dark matter halos have maximum circular velocities in the range 
20 — 35 kms -1 . If the dSphs are embedded within halos this large then the overabundance of satellites 
within the concordance ACDM cosmological model is significantly alleviated, but this still does not 
provide the entire solution. 

Subject headings: cosmology: theory — dark matter — galaxies: halos — halos: evolution — halos: 
structure — methods: iV-body simulations 



1. INTRODUCTION 

The concordance Lambda cold dark matter (ACDM) 
cosmological model for structure formation has been re- 
markably successful in explaining most of the properties of 
the universe on large scales and at various cosmic epochs. 
Recent results from microwave background experiments 
and large redshift surveys have highlighted the ability of 
this model to reproduce observations as diverse as the 
abundance and clustering of galaxies and clusters and the 
statistical properties of the Lya forest within constraints 
set by measurements of the primordial fluctuation spec- 
trum, observations of distant Type la supernovae, and 
gravitational lensing statistics (e.g., Phillips et al. 2001; 
Jaffe et al. 2001; Percival et al. 2001; Hamilton & Tegmark 
2002; Croft et al. 2002; Bahcall et al. 2003; Tonry et al. 
2003; Spergel et al. 2003). However, on non-linear scales 
the ACDM model has been neither convincingly verified 
nor disproved, and several outstanding issues remain un- 
resolved. 

Specifically, the rotation curves of dwarf and low surface 
brightness galaxies are better fitted by shallower dark mat- 
ter density profiles and lower concentrations than those 
predicted by the standard model (e.g., Flores & Primack 
1994; Moore 1994; Burkert 1995; McGaugh & de Blok 
1998; de Blok et al. 2001; de Blok, McGaugh, & Rubin 
2001; de Blok & Bosma 2002; McGaugh, Barker, & de 
Blok 2003; Simon et al. 2003). However, cuspy dark mat- 
ter distributions may be consistent with the rotation curve 
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data once systematic uncertainties are considered, such as 
non-circular motions due to the presence of a bar (Valen- 
zuela & Klypin 2003; Rhee, Klypin, & Valenzuela 2003). 
Fast rotating bars require dark matter densities on galactic 
scales significantly lower than the theoretical predictions 
(Debattista & Sellwood 2000) and fluid dynamical models 
based on observations of barred galaxies support the same 
conclusion (Weiner, Sellwood, & Williams 2001). The in- 
ner Galaxy mass profile is marginally consistent with hav- 
ing dark matter in the central region (Binney & Evans 
2001). On cluster scales, gravitational lensing observa- 
tions are similarly suggestive of a discordance between the 
measured shallow dark matter density inner slopes and the 
predicted by numerical simulations cusps (Sand, Treu, & 
Ellis 2002; Sand et al. 2004). Nevertheless, hydrodynami- 
cal calculations of cluster formation have yet to be carried 
out and compared in detail with the observational data. 

Among the most puzzling discrepancies on small scales 
is the so-called substructure problem. ACDM predicts a 
number of subhalos within the Local Group with maxi- 
mum circular velocities in the range V^a* ~ 10—30 kms -1 , 
which is about 2 orders of magnitude higher than the 
total number of observed satellite galaxies (Kauffmann, 
White, & Guiderdoni 1993; Moore et al. 1999; Klypin 
et al. 1999). A number of possible solutions has been 
proposed to alleviate this problem. One class of solu- 
tions is related to radical modifications of the fundamental 
ACDM paradigm itself, including allowing for a finite dark 
matter particle self-interaction cross section that enhances 
the satellite destruction within galactic halos (Spergel & 
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Steinhardt 2000), reducing the small-scale power with a 
warm dark matter candidate (e.g., Avila- Reese et al. 2001; 
Eke, Navarro, & Steinmetz 2001; Bode, Ostriker, & Turok 

2001) , changing the shape of the primordial power spec- 
trum (Kamionkowski & Liddle 2000; Zentner & Bullock 

2002) , introducing an interaction between dark matter 
particles and photons (Boehm et al. 2002) or resorting to 
the decay of a charged particle to suppress the small-scale 
power spectrum (Sigurdson & Kamionkowski 2003). 

A second class of solutions relies on the inability of low 
mass satellites to form stars, by either supernova feedback, 
photoionization, or rcionization. Indeed, baryon dissipa- 
tion, star formation, and radiative feedback mechanisms 
must have a decisive effect on the properties of the fi- 
nal observed system, and it has been suggested that sup- 
pression of gas collapse and/or cooling by a photoionizing 
background at high redshift can dramatically reduce the 
number of visible satellites with V v i r < 50 kms" 1 , pos- 
sibly reconciling the observations with the model predic- 
tions (Bullock, Kravtsov, & Weinberg 2000; Benson et al. 
2002a, b). Nonetheless, these semi-analytical calculations 
adopt a simple description of the radiation physics and 
disagree with recent numerical simulations showing that 
even very small dwarf galaxies with circular velocities of 
the order of Kir ~ 30 kms -1 or even lower can self-shield 
themselves from the UV background and form some stars 
(Susa & Umemura 2003). Feedback processes may also 
violate the strong correlation between baryonic mass and 
virial velocity, as emphasized by Mayer & Moore (2003). 

Recently, another way of relieving the missing satellites 
problem has been advocated by Stoehr et al. (2002, here- 
after S02). These authors argue that the observed Galac- 
tic satellites are hosted by the most massive substructures 
within a given CDM halo. In this case, feedback pro- 
cesses suppress star formation in all smaller objects, and 
the "luminosity function" of the dwarf spheroidal galaxies 
(dSphs) is consistent with the "mass function" of the sub- 
halos. This represents a somewhat simpler scenario than 
forming stars in a random 10% of the satellites over a larger 
mass range. However, in order to achieve this renormaliza- 
tion of the observed satellite mass/circular velocity func- 
tion, S02 claimed that the observed velocity dispersions 
of the dSphs do not correlate with the maximum circular 
velocities, V max , of their dark matter halos in the way orig- 
inally assumed by Moore et al. (1999) and Klypin et al. 
(1999) when they identified the substructure problem. 

As an illustrative example, consider the case of the dSph 
Fornax, with an observed stellar central velocity disper- 
sion of c7+ <~ 10 kms" 1 (Mateo 1998). The circular ve- 
locity of the surrounding dark matter halo at the loca- 
tion of the dwarf can be of the order of V c = ao* ~ 
15 — 20 kms -1 , depending on the particular assumptions 
made to infer such quantities from observations. For ex- 
ample, an isothermal halo model with a flat circular veloc- 
ity profile, as adopted by Moore et al. (1999), would yield 
V c = V2c± = 14.1 kms -1 . Provided that the circular ve- 
locity profile is still slowly rising in the region where the 
visible dwarf galaxy resides, V c may substantially underes- 
timate V maK . This would clearly allow the dwarf galaxies 
to be embedded in dark halos with very large maximum 
circular velocities in the range of V ma ^ = 50 — 60 kms -1 . 
However, such slowly rising circular velocity profiles seem 



to be in notable disagreement with previous studies sug- 
gesting that the standard ACDM halos corresponding to 
the dSphs are expected to be very concentrated, which 
naturally leads to steeply rising circular velocity profiles 
(Bullock, Kravtsov, & Weinberg 2001). 

S02 suggested that the satellites experience significant 
mass redistribution in their centers as a result of tidal in- 
teractions, leading to shallower inner density profiles and 
smaller effective concentrations than those of comparable 
isolated halos. However, the subhalos in their cosmolog- 
ical simulations contained just a few thousand particles 
and the softening lengths they used were comparable to 
the entire extent of some of the dSphs they wanted to re- 
solve. These authors corroborated their results by compar- 
ing them with the higher resolution TV-body simulations of 
Hayashi et al. (2003, hereafter H03), who employed indi- 
vidual cuspy model satellites disrupting in a static primary 
potential and found the same shallow inner density pro- 
files as a result of tides. Nonetheless, the later models 
have the major shortcoming that they do not constitute 
equilibrium configurations, since they are constructed by 
approximating the exact velocity distribution at any given 
point in space with a Maxwellian. When evolved in isola- 
tion, these models relax rapidly to an inner density slope 
much shallower than the originally intended one, and any 
results obtained should be treated with care (Kazantzidis, 
Magorrian, & Moore 2004). 

In this paper, we investigate how a satellite's internal 
structure responds to tidal interactions, revisiting the so- 
lution to the missing-satellites problem proposed by S02. 
Our study uses A-body cosmological simulations of renor- 
malized systems, in which the subhalos have up to 15 times 
more particles than those used by S02. Note, however, 
that Stoehr et al. (2003) confirmed the results of S02, us- 
ing a simulation with a factor of 9 more particles. Fur- 
thermore, we investigate the tidal disruption of individual 
cuspy substructure halos, employing N = 10 7 particles, 
100 times the mass resolution of H03, to minimize numer- 
ical effects (Moore, Katz, & Lake 1996; Diemand et al. 
2004a). Our primary goal is to examine the change in the 
internal structure of the substructure at scales comparable 
to the actual sizes of the luminous dwarf galaxies as they 
suffer tidal shocks and gravitational stripping. The models 
we adopt for the individual satellite simulations are self- 
consistent realizations and thus ideal for the sensitivity of 
this particular study. As we illustrate below, tides do not 
modify the central density structure of cuspy satellites and 
therefore the missing satellites problem can not be solved 
in this way. 

This paper is organized as follows. In § 2, we discuss 
our cosmological and individual satellite simulations and 
present our results regarding the internal structural evolu- 
tion of the substructure. In § 3, we model the kinematics 
of the Draco and Fornax dSphs on the basis of the findings 
of our simulations, and in § 4 we discuss the implications of 
our results. Finally, in § 5 we summarize our main findings 
and conclusions. 

2. NUMERICAL SIMULATIONS 

All the simulations carried out in this paper were per- 
formed with PKDGRAV (Stadel 2001), a multi-stepping, 
parallel, A-body tree code that uses a spline kernel soften- 



DENSITY PROFILES OF COLD DARK MATTER SUBSTRUCTURE 



3 



ing such that the force is completely Kcplcrian at twice the 
quoted softening lengths. We used an adaptive, kick-drift- 
kick (KDK) leapfrog integrator, and the individual particle 
time-steps A t are chosen according to A t < r](ei/ 'on) 1 / 2 , 
where q is the gravitational softening length of each par- 
ticle, on is the value of the local acceleration, and 77 is 
a parameter that specifies the size of the individual time 
steps and, consequently, the time accuracy of the integra- 
tion. The time integration conserved energy to better than 
0.15% in all cases, which is adequate for the kind of nu- 
merical simulations presented in this paper. The energy 
must be conserved at such a level that the dynamics of 
the regions of interest can be meaningfully probed. We 
have also explicitly checked that our results are not com- 
promised by choices of force-softening, time-stepping, or 
opening angle criteria in the treecode. 

In what follows and unless otherwise explicitly stated, 
we consider as our framework the concordance flat ACDM 
cosmological model with present-day matter and vacuum 
densities il m — 0.3 and £!a = 0.7, respectively, dimen- 
sionless Hubble constant h = 0.7, present-day fluctuation 
amplitude as — 0.9, and index of the primordial power 
spectrum n = 1.0. Note also that in the remainder of 
this paper we use the terms "satellites," "subhalos," and 
"substructure" indistinguishably. 

2.1. Substructure in a Fixed External Potential 

In this section, we investigate the structural evolution 
of CDM substructure, using A-body simulations of satel- 
lites orbiting within the gravitational potential of a static 
primary. The satellites arc modeled using the Navarro, 
Frenk, & White (1996, hereafter NFW) density profile, un- 
der the assumptions of spherical symmetry and isotropic 
velocity dispersion tensors. These models are constructed 
using the procedure described in Kazantzidis ct al. (2004), 
which produces self-consistent equilibria that do not suffer 
from numerical instabilities present in other schemes, such 
as those that approximate the exact velocity distribution 
at any given point in space by a Gaussian. 

The NFW density profile is given by 

= I I ufi TY2 ( r - r ™) ' W 

(r/r s )(l + r/r s y 
where the characteristic inner density p s and scale radius 
r s , are sensitive to the epoch of halo formation and tightly 
correlated with the halo virial parameters, via the concen- 
tration, c = r v i r /r s , and the virial overdensity, A v ; r . Since 
the NFW density profile corresponds to a cumulative mass 
distribution that diverges as r — ► 00, we introduce an ex- 
ponential cutoff for r > r v ; r . The latter sets in at the virial 
radius and turns off the profile on a scale rdecay, which is 
a free parameter and controls the sharpness of the transi- 
tion: 

In order to ensure a smooth transition between equations 
(1) and (2) at r v i r , we require the logarithmic slope there 
to be continuous. This implies 
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The satellites' virial mass is equal to M sat = 1.4 x 
10 10 h~ 1 M Q , corresponding to a circular velocity at the 



virial radius of V\,; r = 35 kms -1 for the adopted ACDM 
model. However, a dark matter halo of a given mass and 
size does not have a unique NFW profile. Indeed, the 
parameter that controls the shape of the density profile 
is the concentration c, and higher values of concentration 
correspond to a larger fraction of the virial mass being con- 
tained in the inner regions of the halo. For a given mass, 
the measured scatter in this parameter reflects mainly 
the different formation epochs (Bullock et al. 2001a; Eke, 
Navarro, & Steinmetz 2001; Wechsler et al. 2002). For the 
adopted ACDM model, the median concentration value 
for an object of this mass scale is c sat = 21, with the 2a 
deviation given by c sat ^11 — 40 (Bullock et al. 2001a). 

Here we present results for two high-resolution simula- 
tions (N = 10 7 ) of subhalos having concentration param- 
eters equal to c sat = 21 (HR1) and c sat = 9 (HR2). We 
consider these two significantly different values of the con- 
centration to single out readily how the evolution of the 
satellites' structure depends on this parameter. We note 
that values as small as c sa t = 9 for the mass scale of our 
satellites are possible within models with a lower fluctua- 
tion amplitude erg. In order to demonstrate the need for 
such a high mass resolution, we run a simulation using a 
lower resolution of N = 5 x 10 5 particles with the lower 
value of the concentration c sat = 9 (LR). We emphasize 
that even the latter run uses a factor of 5 more particles 
than the highest-resolution simulations of H03. 

For the low resolution run, we choose a spline soften- 
ing length of e = 0.51 l/i -1 kpc. This value corresponds 
to an equivalent Plummer softening length equal to that 
of the highest resolution cosmological simulation, GA2, 
of S02. For models HR1 and HR2, we choose a gravita- 
tional softening of e = 0.021/i _1 kpc, comfortably smaller 
than the typical sizes of the dSph satellites of the Milky 
Way which constitute the ultimate target of the present 
study. Numerical and structural parameters of all models 
are summarized in Table 1. For all the runs, we followed 
the time evolution of the density, p(r), and circular veloc- 
ity, V c (r), profiles of the satellites. We explore below how 
these evolve depending on the structural properties of the 
subhalos and on the numerical resolution. 

The orbits of the satellites are influenced only by the 
external, spherically symmetric static tidal field, which is 
represented by the logarithmic halo potential, 

$ = a 2 ln{R 2 + R 2 C ) , (4) 

where a is the one dimensional velocity dispersion and 
R c is the core radius. We use a = 127.3 kms -1 and 
R c = 0.7/i _1 kpc, resulting in a circular velocity profile 
that flattens out already at ~ 3 kpc and reaches an asymp- 
totic value of vo — V~2a ~ 180 kms -1 . This modeling 
gives a total mass within 100 kpc of the center of the fixed 
potential equal to M R<wa kpc = 5.25 x 10 11 h~ 1 M Q . This 
value is at the upper limit of estimates for the mass of the 
Milky Way from satellite motions (Kochanek 1996) and 
dynamics of the Magellanic Clouds (Lin, Jones, & Klemola 
1995), and it is within the range of models for the Milky 
Way proposed in the context of ACDM (Klypin, Zhao, 
& Somervillc 2002). This choice serves our primary aim, 
which is to investigate the change of the satellites' inter- 
nal structure in a regime in which the tidal shocks are very 
strong. Note also that the total mass within the pericenter 
of the orbit is equal to M R<2 5 k pc = 1-3 x 10 11 /i _1 M . 
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Table 1 

Structural and numerical parameters of the satellite models 



Model N c r s V pcak r pcak e sat 

(h- 1 kpc) ( kms- 1 ) (h- 1 kpc) (h- 1 kpc) 
(1) (2) (3) (4) (5) (6) (7) 



LR 


5 x 10 5 


9 


5.4 


41.5 


11.7 


0.511 


HR1 


10 7 


21 


2.3 


51.3 


5.0 


0.021 


HR2 


10 7 


9 


5.4 


41.5 


11.7 


0.021 



Note. — Col. (1): Satellite galaxy model. Col. (2): Number of particles. Col. (3): Concentration parameter. Col. (4): Scale radius. Col. 
(5): Maximum circular velocity. Col. (6): Radius where the circular velocity peaks. Col. (7): Softening length. Note that the virial mass of 
all satellite models is M sat = 1.4 x 10 1() A" 1 M . 



The orbits of the Milky Way dSphs are currently poorly 
constrained observationally. Nevertheless, their current 
distances, which give an indication of the apocenter of 
their orbits, coupled with studies of the orbital proper- 
ties of cosmological halos, can be used to constrain the 
orbital parameters of the satellites. In particular, the 
model satellites are placed on an eccentric orbit with an 
apocenter radius of r apo = 105/i _1 kpc and (r apo /r por ) = 
(6/1), close to the median ratio of apocentric to pericen- 
tric radii found in high resolution cosmological TV-body 
simulations (Ghigna et al. 1998; Tormen, Diaferio, & Syer 
1998). The center of the external potential is always set 
to be (x, y, z) = (0, 0, 0). We begin all our simulations by 
placing the satellites at apocenter, and we integrate the 
orbits forward for 7 Gyr. This timescale corresponds to 
more than three orbital periods (T rb ~ 2.25 Gyr) in the 
chosen orbit and already represents a significant fraction 
of the cosmic time. Our approach neglects the effects of 
the dynamical friction and the response of the primary to 
the presence of the satellite. However, this choice is jus- 
tified given the difference in the mass (almost a factor of 
50) and size of the two systems and the anticipated rapid 
and substantial mass loss owing to tides (see also Taffoni 
et al. 2003). 

The tidal field alters the structure of the satellites through 
a combination of strong tidal shocks occurring at each peri- 
centric passage and gradual mass loss throughout the or- 
bital evolution. In the top panels of Figures 1 and 2, we 
present the evolution of the density (left) and circular ve- 
locity (right) profiles of the bound remnant of the initial 
models HR1 and HR2, respectively. After allowing the 
satellites to pass past the pericenter for the first time, we 
show the profiles at every odd quarter of the orbit between 
apocenter and pericenter. These profiles are plotted from 
the force resolution (2e) outward. The choice of how to 
determine the bound remnants' center is crucial for our 
intended analysis, since quantities such as the spherically 
averaged density and circular velocity profiles are quite 
sensitive to the center's definition. In this study, we iden- 
tify the center adopting the most bound particle method, 



which agrees very well with the center of mass recursively 
calculated using spherical regions of decreasing radius. In- 
deed, we confirmed that both methods give converging re- 
sults by comparing the resulting density profiles for some 
of our stripped satellites. This comparison yielded identi- 
cal profiles for the remnants at all radii between our force 
resolution and the tidal radius for all the cases we tried. 

We determine the mass that remains self-bound as a 
function of time using the following iterative scheme: in 
the rest frame of the most bound particle, we calculate 
the binding energies of all the other particles, using the 
tree-based gravity calculation performed by PKDGRAV, 
and we remove all those with positive binding energy. This 
calculation of binding energies and subsequent removal of 
unbound particles is repeated until no more particles can 
be removed or no bound remnant is found (i.e., all par- 
ticles are removed). In practice, this iterative procedure 
converges rapidly and ensures that the true bound entity 
will be identified. Note that this technique is essentially 
the same one used by most group finding algorithms, such 
as the publicly available SKID (Stadel 2001), which we use 
for the substructure in the cosmological simulations ana- 
lyzed in § 2.2, where removing the background potential 
is more difficult, but it has the advantages of (1) using a 
tree structure for the potential calculation which requires 
of the order of 0(N log N) operations instead of 0(N 2 ), 
where N is the number of particles in the remnant, and 
(2) having a parallel implementation for very large N. In 
this way, we can handle a much larger number of particles 
than would be possible with SKID, and at a fraction of 
the computational burden. 

In the above analysis, we used equal-size logarithmic 
bins. Different number of bins were used, depending on 
the stage of the orbital evolution, ensuring that in each 
case we have a sufficiently large number (> 1000) of parti- 
cles in each bin. This choice of binning simply minimizes 
the noise in the resulting profiles. The central density of 
model HR2 decreases between the initial and final values 
by almost a factor of 3 more than the decrease in model 
HR1 over the same timescales. This is due to the fact that 
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Fig. 1. — Top: Evolution of the density (left) and circular velocity (right) profiles of the bound mass for the high resolution model HR1. 
The number of particles, JV, the concentration parameter, c sa t, and the softening length, e, arc indicated (upper right-hand corner). The 
scale radius of this model is r B = 2.3h~ kpc. The orbital evolution shown corresponds to approximately three orbital periods (T orD ~ 2.25 
Gyr). Profiles are shown at every odd quarter of the orbit between apocenter and pericenter, after allowing the satellites to pass the first 
pcricenter, and are plotted from the force resolution (2e) outward. The lines from top to the bottom, in order of decreasing central density, 
show the profiles at t = (0, 0.67, 1.33, 1.69, 2.36, 2.71) T or b. The density and circular velocity are given in units of M aa t/rf and (GMsat/rs) 1 / 2 , 
respectively. The thick dotted line indicates the expected relation between V max and r max for field halos (see text for details). Bottom: 
Evolution of the central density slope 7 (left) and break radius r D (right) of the fitting formula (eq. [5]) that describes the structure of our 
subhalos. The latter fit parameter is expressed in units of the scale radius, r s , of the initial model. Downward arrows indicate the pericentric 
passages. The substructures maintain the same steep central density slope down to the limit of our force resolution as they get tidally stripped. 
Tides, however, shift inward the break radius of the fitting function. Note the considerable change in rj, after each pericentric passage. 
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Fig. 2. — Same as Figure 1, but for model HR2. The scale radius is r s = 5.4/i — 1 kpc. Compared to model HR1, the decrease in the central 
density and maximum circular velocity, Vmax, is significantly more pronounced, because of the higher binding energy for more concentrated 
models. Tidal interactions truncate this model at smaller physical radii than in model HR1, which is reflected in the evolution of the break 
radius (bottom right). The thick solid line shows the circular velocity profile of the low resolution satellite (LR) at t = 2.71 T or t,. Even 
though this timescale is the same as the one corresponding to the last curve in the high-resolution run, the low resolution satellite has lost 
substantially more mass, and its central slope indicates a shallower inner density profile. Upward arrow indicates the softening used in the 
low resolution simulation. 
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the initial concentrations are quite different. The concen- 
tration parameter of model HR1 is considerably higher 
than that of HR2, and thus the former model is more re- 
silient to tidal heating. The same effect is clearly seen on 
the evolution of the circular velocity profiles. In partic- 
ular, the overall change in the U ma x is significantly more 
pronounced in model HR2 than in model HR1 for the same 
timescales. The thick dotted lines in the top right panels 
of Figures 1 and 2 indicate the expected relation between 
Umax and the distance at which the maximum circular ve- 
locity, r max , occurs for field halos (Colin et al. 2003). The 
latter is a measure of the concentration of the satellites. 
Tidal stripping moves the subhalos to the left of the ex- 
pected relation, so that for a given U max these systems 
have smaller r max and hence higher effective concentra- 
tions than field halos. For example, after <~ 6 Gyr the 
Umax and r max of the c sat = 9 satellite decrease by a fac- 
tor of 3.5 and 8.5, respectively. For the same change of 
Umax, however, one would expect r max to vary by just a 
factor of 3 from the relation for field halos. The higher rel- 
ative densities of heavily stripped halos may account for 
the fact that Draco, being considerably closer to the Milky 
Way, has a mass-to-light ratio higher than that of Fornax, 
despite their similar stellar velocity dispersions. 

The density structure of our subhalos can be described 
through the following simple formula 



p(r) oc r 7 exp( ) 



(5) 



where 7 denotes the central slope of the substructure den- 
sity profile and rb an effective "break" radius describing 
the outer cutoff imposed by the tides. In the bottom pan- 
els of Figures 1 and 2, we show the evolution of these two 
fitting parameters. Note that is expressed in units of 
the scale radius, r s , of the initial models. The robustness 
of the aforementioned binning procedure was verified by 
comparing these results against those when the particles 
were binned in spherical shells with the same number of 
particles in each bin. In this case, we choose as bin center 
the average radius of all particles in the given bin, and 
we assign equal statistical weight to each radial bin. The 
two procedures yielded almost identical values for 7 and rb 
for the entire orbital evolution. Note that in the bottom 
panels of Figures 1 and 2, we present our results adopting 
the latter binning procedure, and the same is also true for 
results shown in Figure 3. 

The bottom left panels of Figures 1 and 2 demonstrate 
that tidal effects do not serve to reduce the central den- 
sity cusp down to the limit of our force resolution. This is 
a fundamental result of this paper that is independent of 
the concentration of the satellites and valid for their en- 
tire orbital evolution. Most of the tidally stripped mass is 
removed from the outer parts of the subhalos, steepening 
the outer density profiles and shifting the break radius rb 
inward. This just reflects the fact that tides truncate the 
satellite models at increasingly smaller physical radii. The 
two bottom right panels illustrate that the lower concen- 
tration satellite (HR2) gets truncated at smaller radii than 
its higher concentration counterpart. In both models, the 
break radii experience the biggest decrease after each peri- 
centric passage, but this decrease is more pronounced in 
model HR2. Between pericenters the break radii decrease 
in a smoother manner. 



Even though this result is also verified in the lower res- 
olution model LR when we compared the density profiles, 
it is interesting to note that the same is not true for the 
circular velocity profiles. Indeed, convergence in the latter 
profile is significantly more problematic, and this is high- 
lighted in the top right panel of Figure 2, where the circular 
velocity profile of the low resolution satellite is shown for 
comparison at t = 2.71 T or b (thick solid line). The down- 
ward arrow indicates the softening used in this simulation. 
The circular velocity profile of a model satellite having a 
factor of 20 fewer particles than its high-resolution coun- 
terpart and evolved on the same orbit is significantly lower 
and has a steeper inner gradient. This is entirely due to 
resolution effects. The fact that the convergence in the 
circular velocity is harder to attain is not surprising, since 
it is a cumulative quantity (see also below, § 2.2). This 
warns against using the circular velocity to compute the 
structural properties of satellites, as done in S02. 

Figures 1 and 2 indicate that the inner satellite regions 
corresponding to sub-kiloparsec scales are much less af- 
fected by the tides than the outer ones. In particular, the 
satellite with c sat = 21 suffers almost no change, even in 
the value of the central density. Subhalo particles that 
have orbital times shorter than the shock duration will be 
only marginally affected by the shock itself. This is known 
as adiabatic correction (Weinberg 1995). Especially for the 
satellite with the highest central density (c sat = 21), the 
adiabatic correction to the impulse approximation is ex- 
tremely important, as it reduces significantly the predicted 
amount of shock heating. We can calculate the variation 
of the kinetic energy of a particle located at some distance 
r from the center of the satellite, AE, due to impulsive 
heating at each pericentric passage in the extended mass 
distribution of a primary isothermal halo (see also Mayer 
et al. 2002; Taffoni et al. 2003). Gnedin, Hernquist, & 
Ostriker (1999) have shown that for orbits as eccentric as 
those considered here, we can approximate the trajectory 
of the satellite with a straight line path and write 



< AE > = 



1 



/ GM > 






UperW 







-2.5 



(6) 

where M and i? ma x are the total mass and radius, respec- 
tively, of the primary halo, which are assumed equal to 
the mass and radius within the apocenter of the satellite's 
orbit, and U pcr is the satellite's velocity at the pericenter 
of the orbit. These parameters are equal to Mo = 8.4 x 
10 11 h^Me, i?max = 105/t- 1 kpc, and U por = 345 kms" 1 . 
The last term in the product is the first-order adiabatic 
correction with w and r being the orbital frequency for a 
particle at radius r and the duration of the tidal shock, 
respectively. 

For a particle at a distance of r = 0.7/i _1 kpc from the 
center of the satellite, AE without the adiabatic correction 
is roughly equal to 20% of its binding energy, -Ebind, after 
three tidal shocks. Instead, A_E/i?bind is already greater 
than 1 at a radius of r = 1.75/i _1 kpc, despite the fact 
that the tidal radius of the satellite is much larger. How- 
ever, once the adiabatic correction is included, the ratio 
A_E/i?bind is reduced by nearly 90%, explaining why the 
satellite is barely affected by the tidal shocks in its inner 
regions. Note that in order to incorporate the effect of the 
adiabatic correction in the simulations, one needs to model 
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Fig. 3. — Density profiles of substructure halos from our individual satellite simulations (solid curves), together with fits to the density 
structure obtained using eq. (5) (circles). The profiles are shown for models HR1 (left) and HR2 (right) and for two different timescales in 
the orbital evolution. In both panels, the density profiles corresponding to the lower curves are vertically shifted downward by 1.5 dex for 
clarity. Curves (a) and (b) correspond to t = 1.69 T or b and ( = 2.71 T or b, respectively. These curves are best described by fitting parameters 
(7> r b/ r s) = (1.02,2.69) and (7,T"b/ r s) = (1.00,2.06), respectively. Curves (c) and (d) correspond to t = 1.33 T or b and t = 2.36T or t, and the 
best fit parameters are (7, r\,/r a ) = (1.02, 0.43) and (7, r^/re) = (0.99, 0.19), respectively. The fits are satisfactory over 4 orders of magnitude 
in density and 2 orders of magnitude in radius and for satellites having concentration parameters that differ by more than a factor of 2. The 
thick dotted line shows the density profile corresponding to the parabolic circular velocity fitting formula proposed by S02 with a = 0.45, the 
median value of their best fits to substructure halos (see text for details). 



the inner regions of the subhalo quite accurately. This 
means resolving down to sub-kiloparsec scales in the first 
place and also avoiding non-equilibrium initial conditions 
or coarse mass resolution. The latter could affect the or- 
bital frequencies of the dark matter particles, which might 
well end up in a region of phase space that is not protected 
by the adiabatic invariance anymore (this problem is rem- 
iniscent of the difficulty that TV-body simulations of low 
resolution have in resolving resonances in stellar systems; 
see Weinberg & Katz (2002)). 

Figure 3 shows the density profiles for two different 
timescales in the orbital evolution of models HR1 and 
HR2 (solid curves). Fits to the density structure of the 
substructure halos using equation (5) and obtained by x 2 
minimization are also shown, as circles. Model HR1 (left) 
is plotted at t = 1.69 T orb and 2.71 T" rb, which correspond 
to ~ 0.4 Gyr after the satellite has concluded its second 
and third pericentric passages, respectively. Model HR2 
(right) is plotted at t = 1.33 T orb and 2.36 T orb , which cor- 
respond to ~ 0.4 Gyr before the conclusion of the satel- 
lite's second and third pericentric passages, respectively. 
Note that for clarity we offset the lower curves in both 
panels by 1.5 dex. Clearly, the structure of our stripped 
satellites is reasonably well reproduced by the proposed 
fitting formula. We should note that this formula works 
well in the case of the lower concentration satellite for the 
entire orbital evolution. On the other hand, the case of the 
c sat = 21 satellite is somehow slightly different. The fitting 
function works better for the last stages of the orbital evo- 
lution, at which the satellite has suffered significant mass 
loss. This is just reflecting the resilience to tidal stripping 
of the high concentration satellite during the early stages 
of the evolution and the fact that the mass loss is more 
gradual compared to the low concentration model, even 
from the outer regions. Finally, for comparison we show 



the density profile resulting from the fitting formula pro- 
posed by S02 (thick dotted line) to describe the structure 
of their subhalos (see their eq. [1]), plotted down to the 
limit of our force resolution. The values for r max and V max 
in their formula are taken directly from the particular sub- 
halo. In addition, we adopt a value for the parameter a 
in their equation that is equal to the median of their best 
fits, a — 0.45. Compared to the density profiles describ- 
ing the structure of our satellites, the latter profile has a 
substantially shallower slope on scales comparable to the 
sizes of the dSphs. Interestingly, the density reaches even 
negative values at a finite radius, demonstrating that its 
use should be avoided. 

2.2. Substructure in Hierarchical Cosmological 
Simulations 

In this section, we present results from a set of high reso- 
lution ACDM cluster simulations (Diemand et al. 2004b). 
The initial conditions are generated with the GRAFIC2 
package (Bertschinger 2001). We begin with a 300 3 par- 
ticle cubic grid with a comoving cube size of 300 Mpc 
(particle mass m p = 2.6 x lO lo /i _1 M ). At z = we 
identify several clusters for resimulating at higher resolu- 
tion, using refinement factors of 6, 9, and 12 in length 
(216, 729, and 1728 in mass), so that the mass resolu- 
tion is to p = 1.5 x 10 7 /i _1 M Q in the highest resolution 
run. We label these three runs as R6, R9 and R12, af- 
ter their refinement factors. At z = the refined cluster 
contains 1.8 x 10 6 particles within the virial radius in run 
R6, 6 x 10 6 in R9, and 1.4 x 10 7 in R12. Note that we 
define the virial radius of the clusters, i? v ; r , to be the dis- 
tance from the center at which the mean enclosed density 
is 178f2^ 45 w 103.5 times the critical value p cr i t (Eke, Cole, 
& Frenk 1996). The softening length is comoving from the 
start of the simulation (z ~ 40) to z = 9. From z — 9 until 
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Fig. 4. — Top: Density (left) and circular velocity (right) profiles of two subhalos simulated at two different resolutions. Circular velocities 
are expressed in units of the maximum circular velocity of the parent halo. In simulation R9 (solid lines), the subhalos contain 3.4 X 10 4 
and 1.4 X 10 4 particles, while in R6 (dashed lines) they contain a factor of 3.375 fewer particles. The convergence in the circular velocity 
profiles between the two resolutions is significantly slower on account of the cumulative nature of this quantity. This demonstrates that it 
is erroneous to use the circular velocity to compute the structural properties of substructure. Bottom left: Density profiles of two heavily 
stripped substructure halos in simulation R12 before entering the primary halo (dashed lines) and at present (solid lines). Subhalo (i) is 
shown at z = 0.80 and just before the second pericentric passage. Subhalo (ii) is shown at z = 0.94 and just before the fourth pericentric 
passage and is offset by 1 dex to avoid overlap. The numbers near the thick solid lines indicate the power slope of those lines. Both subhalos 
have of the order of N = 2 X 10 4 particles before entering the host. The cosmological subhalos maintain their steep inner density slope even 
after several pericentric passages. Bottom right: Density profiles of the five most massive halos in simulation R12 (solid lines). The best-fit 
NFW profiles are also plotted (dashed lines). To avoid overlap the lower four density profiles are vertically shifted by 1, 2, 3, and 4 dex from 
the top to the bottom. These five halos are resolved with more than 2 X 10 5 particles. In all panels the vertical dotted lines show the adopted 
force resolution. 
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the present we use a physical softening length of 2.1, 1.4, 
and 1.05 x 10 _3 i? v ir for R6,R9, and R12 respectively. 

In the top panels of Figure 4, we show the density {left) 
and circular velocity (right) profiles of two distinct sub- 
structure halos simulated at two different resolutions. In 
simulation R9 (solid lines), the two subhalos contain N = 
3.4 x 10 4 and N = 1.4 x 10 4 particles, respectively, whereas 
in simulation R6 (dashed lines) the same subhalos con- 
tain a factor of 3.375 fewer particles. The vertical dotted 
lines show the adopted force resolution. Densities are ex- 
pressed in units of the present critical density for closure, 
Pcrit = 3i?o /87rG, and the virial radius of the parent halo 
is equal to i? v ir = 1.225ft- -1 Mpc. The profiles of the cos- 
mological satellites are to be trusted only up to the res- 
olution limit set by two body relaxation (Diemand ct al. 
2004a). We note that the convergence of the circular ve- 
locity profiles is much weaker than that of the density pro- 
files at different resolutions, in agreement with the results 
presented in § 2.1. As the circular velocity is a cumula- 
tive quantity, numerical effects show up at comparatively 
larger radii, where the density profiles have already con- 
verged, biasing the circular velocity profiles. Therefore, 
one should be especially cautious when using the circular 
velocity profiles obtained at one particular resolution to 
provide fitting functions to systems whose typical scales 
fall below such resolution, as done in S02. 

In the bottom panels of Figure 4, we present results re- 
garding substructure halos taken from our highest resolu- 
tion simulation (R12). The density profiles of two heavily 
stripped subhalos before entering the primary halo (dashed 
lines) and at present (solid lines) are shown in the bottom 
left panel. These satellites have of the order of N — 2 x 10 4 
particles before entering the host. The subhalo denoted 
"(i)" is shown at z; n = 0.80 and at present just before the 
second pericentric passage (r m i n = 0.21i? v i r ). The second 
subhalo, denoted "(ii)," is offset by 1 dex to avoid overlap 
and is shown at z ln = 0.94 and at present just before the 
fourth pericentric passage (r m - m — 0.15i? v ir)- The thick 
solid lines indicate two critical values for the inner slope, 
p(r) oc r -1 and p(r) oc r -1 ' 5 . This plot illustrates that the 
substructures maintain the steep inner density slope that 
they had before entering the primary halo, even after sev- 
eral pericentric passages. These heavily stripped satellites 
have outer profiles consistent with an exponential cutoff. 
Even though these cosmological simulations are state-of- 
the art by the current standards, the resolution in these 
stripped satellites is still too low for us to draw robust con- 
clusions regarding the evolution of their inner structure. 

Finally, the bottom right panel of Figure 4 shows the 
density profiles of the five most massive substructure ha- 
los in this simulation (solid lines) together with the best fit 
NFW profiles (dashed lines). To avoid overlap, the lower 
four density profiles are each vertically shifted by 1 dex 
from each other. These five halos are resolved with more 
than N = 2 x 10 5 particles, and their density profiles are 
described reasonably well by the NFW profile. These halos 
entered the host system quite late and have typically made 
only one orbit; however, these most massive halos repre- 
sent the ones that would be associated with the dSphs, 
according to S02. We can plausibly trust their central 
structure to 0.002i? v j r , which corresponds to about 500 pc 
when scaled to the Galaxy. At this radius the density pro- 



files are still cuspy, with central slopes between —1 and 
-1.5. 

We end by emphasizing that the behavior of the satel- 
lites in the time-dependent cosmological tidal field is not 
obviously similar to that our individual subhalos exhibit 
when evolved in a static host potential. The latter are 
spherical systems with isotropic velocity dispersion ten- 
sors, whereas the cosmological satellites presented here 
are in general triaxial, anisotropic systems, and they suffer 
additional artificial heating from background particles and 
physical heating from encounters with other substructures. 
It is remarkable, therefore, that our findings regarding the 
way that tidal interactions affect the central density cusps 
converge in these two radically different regimes: the inner 
density profiles remain cuspy to the resolution limit of the 
simulations. 

3. DWARF SPHEROIDAL KINEMATICS IN CDM SUBHALOS 

Cosmological numerical simulations confirm that iso- 
lated halos in the mass range 10 7 — 10 10 M Q have central 
cusps not shallower than p(r) oc r , on scales compara- 
ble to the ones probed by the stars in the dSph satellites 
(Moore et al. 2001; Colin ct al. 2003). Our results in § 2 
indicate that an initial cuspy satellite will retain its cusp in 
the presence of a tidal field. Therefore, it is an interesting 
exercise to attempt to reproduce the observed kinematics 
of the dSphs assuming cuspy dark matter distributions, 
and to constrain the maximum mass or circular velocity 
of their host halos. The advantage of our approach lies in 
the fact that we have sufficient numerical resolution, and 
therefore there is no need of extrapolating the inner den- 
sity slopes to scales smaller than the force resolution, as 
done by S02. 

The dynamics of a spherically symmetric stellar distri- 
bution embedded in a spherical dark matter potential can 
be described by the lowest order Jeans equation, assuming 
that the two components are in equilibrium and that there 
are no net streaming motions (e.g., rotation): 

— (pa I ) + —pa I +p— = (7) 

(Binney & Tremaine 1987), where p(r) and cr r (r) are the 
density profile and the radial velocity dispersion of the 
tracer stellar population, respectively, <&(r) is the under- 
lying dark matter gravitational potential, and j3 describes 
the velocity anisotropy of the stars. The solution of the 
Jeans equation (7) with the boundary condition p(r) — > 
at r — > r t for /3=const reads 

pn^ = r- 2 ' 3 P r 2 Pp^- dr , (8) 
J r dr 

where r t is the tidal radius of the stellar component. 

However, the quantity that can be measured and com- 
pared with observations is the line-of-sight velocity dis- 
persion, cjp, at projected distance R from the center of the 
dwarf, which is given by 

(Binney & Mamon 1982; Binney & Tremaine 1987), where 
I(R) is the surface distribution of the tracer population: 
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Fig. 5. — Observed line-of-sight velocity dispersions for the dSph Draco (triangles), compared with those predicted for stellar systems 
embedded in dark matter halos with structure similar to the ones of our simulated subhalos. Each pair of panels refers to a single value of 
Vmax, while the various lines correspond to different values of Rmax in the same range as in S02 (see text for details). The three plots on the 
left show results assuming a tidal radius for Draco equal to the optical radius measured by Odenkirchen et al. (2001) and also adopted by S02. 
The panels on the right show results for the same range of parameters, but for tidal radii 3 times as large. The filled points correspond to the 
results obtained assuming the parabolic fits to the circular velocity profiles adopted by S02. From top to bottom, the results for Vmax = 54, 
35, and 25 kms -1 are shown. Globally, the observed velocity dispersion profile is better reproduced by subhalos with V ma x = 25 kms -1 and 
assuming a significantly larger than the nominal tidal radius. 



We solve equation (9) using the customary King profile 
for the stellar population of the dSph galaxies, 



p*(r) oc 



1 



icos- 1 (z)-(l-z 2 ) 1 / 2 



(11) 



(King 1962), where z = [1 + (r/r c ) 2 ]/[l + (r t /r c ) 2 ] and r c 
denotes the core radius of the profile. For Fornax and 
Draco, we use the parameters from Mateo (1998) and 
Odenkirchen et al. (2001), respectively, similarly as in S02. 
Note that the normalization in equation (11) is not im- 
portant for the purpose of our analysis, since it cancels 
out. Motivated by the structure of our simulated subha- 
los from the cosmological simulations, we assume that the 
dark matter halos associated with the dwarf galaxies fol- 
low the NFW density profile. This provides a lower limit 
to the recovered velocity dispersions, since there is evi- 
dence in favor of steeper central density cusps on these 
scales (e.g., Moore et al. 2001). Our results are not sen- 
sitive to the presence of the exponential cutoff, since the 
fitting formula of equation (5) differs appreciably from the 
NFW density profile only near the break radius ?v This is 
larger than the tidal radii of the stellar components of the 
dSphs, where the integral in equation (9) has to be trun- 
cated anyway. In addition, we checked that any eventual 
differences in the resulting projected stellar velocity dis- 
persions were negligible even when we considered models 
with stellar tidal radii significantly larger than the nominal 
values (see below). 

We consider halos having both a maximum circular ve- 
locity, Umax, and a corresponding radius, i? max = R(V max ), 



within the range considered by S02. In practice, for each 
value of U m ax there is only a limited range of halo con- 
centrations such that i? max is consistent with the values 
measured by S02 for their subhalos. In Figures 5 and 6, 
we show the results of the comparison for Draco and For- 
nax, respectively, assuming isotropic stellar orbits. The 
observed velocity dispersion profiles are reproduced from 
Mateo (1998) for Fornax and from Kleyna et al. (2002) for 
Draco. In some of these panels, where the relative informa- 
tion is available, we overplot the results obtained adopting 
the parabolic fitting function of the form proposed by S02 
for their subhalos (filled points). We remind the reader 
that the density profiles corresponding to this parabolic 
function have substantially shallower inner slopes at scales 
comparable to the sizes of the dSphs than the profiles de- 
scribing the density structure of our satellites. Each pair 
of panels refers to a single value of U ma x- The three panels 
on the left show results assuming a tidal radius equal to 
that adopted by S02, while the panels on the right assume 
tidal radii 3 times as large. 

It is immediately apparent that subhalos with Umax ~ 
50 km s _1 would yield velocity dispersions significantly above 
the observed values for both Fornax and Draco, once the 
cuspy density profiles consistent with the structure of our 
satellites were used. The disrepancy becomes even worse if 
we allow the stellar component to have a tidal radius larger 
than the nominal value estimated by fitting the star counts 
to a King profile. Recent observational studies of Draco 
show no evidence of tidal tails, lending support to the view 
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Fig. 6. — Same as Figure 5, but for the dSph Fornax (squares). The three plots on the left show results assuming a value for the tidal 
radius taken from the review of Mateo (1998). From top to bottom, the results for Vmax = 46, 31, and 25 kms" 1 are shown. Similarly to the 
case of Draco, the stellar velocity dispersion is better reproduced by subhalos with Vmax = 25 kms -1 and assuming a tidal radius 3 times 
larger than the nominal value (right). 
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Fig. 7. — Kinematics of Fornax (top) and Draco (bottom), compared with those expected for stellar systems embedded within NFW 
subhalos with Vmax = 25 kms" 1 . Different values for the anisotropy parameter /3 in the velocity distribution of the stellar component have 
been considered. The solid lines correspond to isotropic models, whereas the dotted and dashed lines to radially and tangentially anisotropic 
models, respectively. The plots on the left show results for values of the tidal radii of the stellar component equivalent to those used by S02, 
while the plots on the right correspond to tidal radii 3 times as large. Models with mildly tangentially anisotropic distribution of stellar orbits 
better reproduce the observed velocity dispersion profiles of both Draco and Fornax. Radially anisotropic models overestimate the central 
stellar velocity dispersion and make the curves decrease more steeply with distance, contrary to the trend in the observational data. 
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Fig. 8. — Best fit models (solid lines) and 3cr intervals (dashed lines) for the kinematics of Draco, assuming isotropic (top) and tangentially 
anisotropic (bottom) models for the stellar distribution. The stars are embedded in the potential wells of our stripped satellites and follow a 
King profile. The panels on the left show results adopting a tidal radius equal to the nominal value, whereas the panels on the right show 
results for tidal radii 3 times as large (see text for details). Keeping the concentration parameter fixed (c = 21), wc fit by \ 2 minimization 



the maximum circular velocity V ma x. The best fit value of V max corresponds to v be 



whereas 



and v n 



bracket the 3<r intervals. The 



quality of the fits in terms of the reduced \ is also indicated for the different cases. Subhalos with Vmax in the range 20 - 
the best fit to the data, while Vmax > 40 kms are 3<r or more away from the best fit. 



35 kms provide 



that the optical radius is only a lower limit to the true 
physical boundary of the system (Odcnkirchen et al. 2001; 
Kleyna et al. 2002; Piatek et al. 2002; Klessen, Grebel, & 
Harbeck 2003), and this might also be true for the major- 
ity of the dSphs. By considering subhalos with U max <~ 
30 kms -1 , we match satisfactorily only the velocity dis- 
persion data points at the outermost radii, which have the 
largest error bars. Subhalos with U max = 25 kms -1 repro- 
duce the observed stellar velocity dispersions reasonably 
well for the entire range of parameters considered. This 
conclusion is in disagreement with that reached by S02, 
simply because of their use of a parabolic fitting function, 
and suggests that there is no room for entirely solving 
the substructure problem by simply changing the way the 
velocities are mapped. Globally, the observed velocity dis- 
persion profile is better reproduced by assuming a larger 
tidal radius, even though the predicted profiles always tend 
to be slightly flatter than the data. 

The situation can be improved significantly by adopting 
a weakly tangentially anisotropic velocity distribution for 
the stars, as shown in Figure 7. Indeed, for /3 = — 1 (dashed 
lines), subhalos with U max = 25 kms -1 arc nicely consis- 
tent with most of the data points for both Draco and For- 
nax. Note that a mildly tangentially anisotropic distribu- 
tion of stellar orbits was required by Kleyna et al. (2002) in 
their best fit models to the kinematics of Draco. Moreover, 
Lokas (2002) finds that the steeper the dark halo inner den- 
sity cusp, the more tangential a stellar velocity distribution 
is required to fit the data for both Draco and Fornax. In- 
terestingly, modeling based on high-resolution Keck spec- 



troscopy of six Virgo dwarf elliptical galaxies having an 
average line-of-sight velocity dispersion in the range a ~ 
25 — 50 km s -1 also yields for most of them best fit models 
to their kinematics that are consistent with mildly tangen- 
tially anisotropic orbits (Geha, Guhathakurta, & van der 
Marel 2002). Radially anisotropic models with (i = 0.5 
(dotted lines) overestimate the central stellar velocity dis- 
persion by a factor of ~ 2 and lead to steeply declining 
velocity dispersion profiles, contrary to the trend in the 
observational data. As in the case of the isotropic orbits, 
if we assume a tidal radius 3 times larger (right) the match 
to the stellar velocity profiles becomes better. 

In Figures 8 and 9, we present the best fit models (solid 
lines), together with the 3cr intervals (dashed lines) for 
Draco and Fornax respectively, assuming isotropic (top) 
and tangentially anisotropic (bottom) models. In each 
case, we fix the value of the concentration parameter and 
fit, by x 2 minimization the maximum circular velocity 
Knax with the two different values for the tidal radius of 
the stellar component. Note that for Draco the best fit 
is calculated taking into account that the error bars are 
asymmetric. The fixed concentration parameter we used 
(c = 21) is consistent with nearly all the profiles of the 
subhalos in S02. Higher concentration models would im- 
ply lower values of U max , and therefore our choice is con- 
servative. This more detailed analysis shows that subhalos 
with U max in the range 20 — 35 kms -1 provide the best fit 
to the data and also that the high values favoured by S02, 
U max > 40 kms -1 , are 3cr or more away from the best fit. 
Determining the absolute best fit to the data requires fit- 
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Fig. 9. — Same as Figure 8, but for the dSph Fornax. Similarly to the case of Draco, subhalos with Vmax in the range 20 — 30 kms 
provide the best fit to the data (solid lines), while Vmax > 35 kms -1 are 3<r (dashed lines) or more away from the best fit. 



ting more parameters, which is clearly outside the scope of 
the present paper (however, see Lokas 2002). In general, 
the degree of anisotropy can be a function of radius instead 
of being simply a constant parameter. In models in which 
the dSphs are formed by the tidal stirring and transforma- 
tion of disk-like systems similar to present-day dwarf irreg- 
ular galaxies (Mayer et al. 2001a, b), the remnant is a tri- 
axial system whose velocity anisotropy depends on the ra- 
dius. The remnants produced in those simulations can be 
fitted by cither an exponential or a King profile, similar to 
the observed dSphs. Although there is considerable scat- 
ter in their structural properties (Mayer et al. 2001b), in 
general the stellar orbits are tangentially anisotropic in the 
central regions (at distances equivalent to the core radii of 
Draco and Fornax) and become nearly isotropic or slightly 
radially anisotropic in the outer regions, where they are 
partially flattened by rotation. The interesting point is 
that values of (3 ranging from —0.5 to —1.3 are typical in 
the region where most of their bound mass is. This result 
is confirmed by recent gas-dynamical simulations with ra- 
diative cooling and heating and star formation (L. Mayer 
et al. 2004, in preparation; Mayer & Wadsley 2003) and 
lends support to the results shown in Figure 7 for (i = — 1 . 
Note that Zentncr & Bullock (2003) suggested that a mild 
radial anisotropy for the dSphs equal to (3 = 0.15 should be 
expected because the latter value is typical of the central 
regions of simulated CDM halos, where the dwarfs actu- 
ally reside. However, our results highlight that this is not 
necessarily the case. The tidally stirred stellar components 
develop a tangential anisotropy even though they evolve 
inside an initially isotropic dark matter halo because their 
evolution is governed by non-axisymmetric, bar/buckling 
instabilities (Mayer et al. 2001b). 

We end by emphasizing that if the density profiles of the 



subhalos models are initially steeper than an NFW profile 
by either having intrinsically steeper inner slopes (as might 
be possible for some of the satellites in our cosmological 
runs; see § 2.2) or becoming steeper as a result of the adi- 
abatic contraction following the collapse of baryons (Blu- 
mcnthal ct al. 1986), circular velocities even lower than 
the above values would be required to acceptably repro- 
duce the observed data. 

4. DISCUSSION 

The results presented in the previous sections indicate 
that the conversion between ov and K max initially adopted 
by Moore et al. (1999) and Klypin et al. (1999) is rea- 
sonable. Low values of Vmax are required if the dSphs 
are embedded in CDM halos. There is a second, phe- 
nomenological reason why the observed satellites should 
not reside in halos as massive as the ones postulated by 
S02. In this case, their baryonic components would lie 
within a very small region corresponding to only 1% of the 
virial radius of the object, posing a problem for structure 
formation scenarios. Indeed, in the prevailing galaxy for- 
mation paradigm (e.g., Fall & Efstathiou 1980; Mo, Mao, 
& White 1998), the baryons condense into a rotationally 
supported structure whose size is determined by the initial 
spin parameter A and concentration c of the dark matter 
halo, and the fraction of mass and angular momentum of 
baryons relative to that of the halo. Assuming that the 
specific angular momentum of baryons is conserved dur- 
ing their infall and that the halo and baryons start with 
the same specific angular momentum (Mo, Mao, & White 
1998), one can easily show that for a halo as massive as 
V v ir = 50 kms -1 , a spin parameter equal to A < 0.01 
would have been needed for the baryons to collapse to the 
observed size of Draco. However, such a small value for 
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the spin occurs in less than 1% of the cosmological halos 
(Warren et al. 1992; Lemson & Kauffmann 1999; Gardner 
2001). 

The distribution of spin parameters has been measured 
so far only for isolated halos and at mass scales larger 
than those of dwarf galaxies (e.g., Bullock et al. 2001b; 
Gardner 2001). Although a new systematic analysis over- 
coming the above limitations will be needed in the future, 
we made the first step forward in this direction and we 
measured the spin parameter A for some of our subhalos. 
In particular, we selected ten isolated halos from the high 
resolution region of run R6 and six subhalos within the 
virial radius of the cluster in both runs R6 and R12. All 
halos are selected by mass (M v ; r ~ 6.3 x 10 10 ft~ 1 M Q ). 
The spin parameters are measured at a radius equal to 
0.5r v j r for the isolated halos, where r v ; r denotes the halo 
virial radius, and at 0.7 r t for the subhalos, where r t de- 
notes the subhalo tidal radius. The A-values we found for 
both isolated halos and subhalos are between 0.024 and 
0.12 and follow a log-normal distribution with A = 0.039 
and c \ = 0.55 in agreement with previous studies using a 
much larger halo sample (e.g., Bullock et al. 2001b). This 
finding confirms our phcnomcnological argument that the 
dSph satellites could not be embedded in very massive ha- 
los. Finally, it is interesting to note that the comparison 
of the spin parameters of the same subhalos in runs R6 
and R12 shows that they are independent of resolution. 
Removal of baryons from an initially more extended sys- 
tem by tidal stripping, by supernova feedback, or even by 
photoevaporation from the UV background cannot be in- 
voked to reduce the size of the baryonic component of the 
dwarfs at a later stage, because the halo would have been 
massive enough to suppress all these mechanisms (Ben- 
son et al. 2002b). One then would have to rely on some 
catastrophic loss of angular momentum during baryonic 
collapse, but there is no obvious reason why this should 
have happened. 

Recently, Kravtsov, Gnedin, & Klypin (2004) proposed 
a qualitatively different solution to the missing satellites 
problem. Using high resolution cosmological simulations 
of the formation of a Milky Way sized halo in a ACDM 
universe, these authors argued that the luminous dSphs in 
the Local Group can be identified with halos that had con- 
siderably higher masses and circular velocities when they 
formed at high redshift. Their model provide a convinc- 
ing explanation as to why the dSphs could retain their gas 
and form stars after reionization and indicates that they 
are embedded in dark matter halos that suffered dramatic 
mass loss due to tidal stripping. We believe that the major 
difference between the subhalos' density structure found in 
our study and that claimed by S02 is due to resolution ef- 
fects affecting the determination of the circular velocity. 
On the other hand, the reduced concentration that H03 
find for their simulated satellites may be an artifact of 
non-equilibrium initial conditions combined with numeri- 
cal resolution, rather than reflecting the influence of tidal 
shocks on the inner regions of the satellites. The evolution 
of the internal structure is expected to be dramatically dif- 
ferent for models that are not self-consistent compared to 
models for which one considers carefully the exact dynam- 
ics. Regarding this issue, Kazantzidis et al. (2004) found 
that NFW satellites initialized using the Maxwellian ap- 



proximation experience artificially accelerated mass loss 
and can completely disrupt in an external tidal field. The 
same satellites, however, survive in identical experiments 
once the exact self-consistent velocity distribution of the 
model is taken into account. 

Of course, one cannot exclude that some dynamical mech- 
anism other than tidal interactions might have been re- 
sponsible for lowering the concentration of the subhalos. 
For example, if bars are effective in redistributing the an- 
gular momentum between the baryons and the halo (Wein- 
berg & Katz 2002), then the dSphs could have been strongly 
affected by this mechanism if they originate from the mor- 
phological evolution of tidally stirred dwarf irregular galax- 
ies. Indeed, tidal stirring always incorporates a strong bar 
instability phase. However, recent work suggests that the 
amount of angular momentum that a bar transfers to the 
halo is not enough to turn the inner cusp into a core (Sell- 
wood 2003). 

The results presented in this paper argue that there is no 
simple solution to the substructure problem within CDM 
models. The kinematics of the dSphs are indeed a rea- 
sonable tracer of their dark halo potential wells. However, 
their kinematics can also be accurately reproduced with 
dark halo density profiles that are nearly flat in the center 
(see Figures 5 and 6). The latter fact illustrates the need 
to use better quality data and perhaps higher order mo- 
ments of the velocity distribution to break this degeneracy 
and constrain the central cusp slopes (Lokas 2002). Our 
findings also have important implications for indirect dark 
matter detection experiments attempting to observe the 
resulting 7-ray flux from neutralino annihilation. The ex- 
pected flux depends sensitively on the line-of-sight integral 
of the square of the mass density. This integral changes 
by more than one order of magnitude between the NFW 
and the much shallower density profiles found by S02, as 
pointed out by the same authors in a more recent paper 
(Stoehr et al. 2003). With such shallow density profiles, 
Air-Shower-Cerenkov telescopes like Veritas and Magic 
might barely detect the emission from the Galactic center 
and would completely miss that coming from the dSphs. 
However, for dark matter particle candidates with proper- 
ties consistent with some of the minimal supersymmctric 
models (Bergstrom et al. 1998) and with the persistently 
steep inner density profiles found in this paper, the Galac- 
tic center would be easily detected and the dSph satel- 
lites still represent potentially valuable targets (Calcaneo- 
Roldan & Moore 2000). 

5. SUMMARY 

We have investigated the structural evolution of sub- 
structure, using a set of high-resolution cosmological N- 
body simulations, coupled with simulations of the tidal 
stripping of individual satellites orbiting within a static 
host potential and employing up to N = 10 7 particles. 
Our main results and conclusions may be summarized as 
follows: 

1. Cuspy satellite halos on eccentric orbits do not ex- 
perience significant mass redistribution in their cen- 
ters, and they retain the same steep inner density 
slope even after several strong pericentric tidal shocks. 
This result is valid for our high-resolution cosmologi- 
cal simulations, in which the satellites evolve within 
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a time-dependent cosmological tidal field, and for 
our simulations of individual subhalos orbiting within 
a massive host system modeled fixed potential. 

2. Convergence in the circular velocity profiles occurs 
much more slowly than that of the density profiles. 
This is due to the fact that low central resolution 
propagates to larger radii on account of the cumula- 
tive nature of the former quantity. This discourages 
any attempt to use circular velocities of simulated 
satellites to derive their structural properties. 

3. The density structure of tidally stripped cuspy ha- 
los can be well approximated by a simple fitting 
function (eq. [5]), comprising an unmodified power- 
law central slope and an exponential cutoff, that 
describes the satellite's boundary imposed by the 
tides. 

4. The predicted kinematics of the dSph galaxies For- 
nax and Draco, assuming that they are embedded in 
the potential wells of our tidally stripped satellites, 
were compared to the observed stellar velocity dis- 
persion profiles. Adopting isotropic and tangcntially 
anisotropic models for the stellar component, we 
find that dark matter halos with maximum circular 
velocities in the range V max <~ 20 — 35 kms -1 fit the 
data better, whereas models with V max & 40 kms -1 
are at least at the 3<r level from the best fit values. 
If the initial central density slopes of the simulated 
subhalos were steeper than the ones adopted here, 
even lower values of Vmax would be required to fit the 
data. Tidal interactions do not provide the mecha- 
nism for embedding the luminous dwarf galaxies of 
the Milky Way within the most massive subhalos 
in a ACDM universe, and, therefore, for reconciling 
the overabundance of Galactic satellites with CDM 
predictions. 

5. Models with a mild tangential anisotropy of stellar 
orbits reproduce better the shape of the observed ve- 
locity dispersion profiles of both Fornax and Draco, 
whereas radially anisotropic models overestimate the 
central stellar velocity dispersion by a factor of ~ 2. 
Such a tangential anisotropy is expected in scenarios 
in which the dSphs result from the tidal stirring of 
systems similar to present-day dwarf irregular galax- 
ies. 
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